home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / public / pbmplus / pgm / fitstopgm.c < prev    next >
C/C++ Source or Header  |  1994-08-01  |  6KB  |  242 lines
  1. /* fitstopgm.c - read a FITS file and produce a portable graymap
  2. **
  3. ** Copyright (C) 1989 by Jef Poskanzer.
  4. **
  5. ** Permission to use, copy, modify, and distribute this software and its
  6. ** documentation for any purpose and without fee is hereby granted, provided
  7. ** that the above copyright notice appear in all copies and that both that
  8. ** copyright notice and this permission notice appear in supporting
  9. ** documentation.  This software is provided "as is" without express or
  10. ** implied warranty.
  11. */
  12.  
  13. #include "pgm.h"
  14.  
  15. struct FITS_Header {
  16.     int simple;        /* basic format or not */
  17.     int bitpix;        /* number of bits per pixel */
  18.     int naxis;        /* number of axes */
  19.     int naxis1;        /* number of points on axis 1 */
  20.     int naxis2;        /* number of points on axis 2 */
  21.     int naxis3;        /* number of points on axis 3 */
  22.     double datamin;    /* min # */
  23.     double datamax;    /* max # */
  24.     double bzer;    /* ??? */
  25.     double bscale;    /* ??? */
  26.     };
  27.  
  28. static void read_fits_header ARGS(( FILE* fp, struct FITS_Header* hP ));
  29. static void read_card ARGS(( FILE* fp, char* buf ));
  30.  
  31. void
  32. main( argc, argv )
  33. int argc;
  34. char* argv[];
  35.     {
  36.     FILE* ifp;
  37.     gray* grayrow;
  38.     register gray* gP;
  39.     int argn, imagenum, image, row;
  40.     register int col;
  41.     gray maxval;
  42.     double fmaxval, scale;
  43.     int rows, cols, images;
  44.     struct FITS_Header h;
  45.     char* usage = "[-image N] [FITSfile]";
  46.  
  47.     pgm_init( &argc, argv );
  48.  
  49.     argn = 1;
  50.     imagenum = 1;
  51.  
  52.     if ( argn < argc && argv[argn][0] == '-' && argv[argn][1] != '\0' )
  53.     {
  54.     if ( pm_keymatch( argv[argn], "-image", 2 ) )
  55.         {
  56.         ++argn;
  57.         if ( argn == argc || sscanf( argv[argn], "%d", &imagenum ) != 1 )
  58.         pm_usage( usage );
  59.         }
  60.     else
  61.         pm_usage( usage );
  62.     ++argn;
  63.     }
  64.  
  65.     if ( argn < argc )
  66.     {
  67.     ifp = pm_openr( argv[argn] );
  68.     ++argn;
  69.     }
  70.     else
  71.     ifp = stdin;
  72.  
  73.     if ( argn != argc )
  74.     pm_usage( usage );
  75.  
  76.     read_fits_header( ifp, &h );
  77.  
  78.     if ( ! h.simple )
  79.     pm_error( "FITS file is not in simple format, can't read" );
  80.     switch ( h.bitpix )
  81.     {
  82.     case 8:
  83.     fmaxval = 255.0;
  84.     break;
  85.  
  86.     case 16:
  87.     fmaxval = 65535.0;
  88.     break;
  89.  
  90.     case 32:
  91.     fmaxval = 4294967295.0;
  92.     break;
  93.  
  94.     default:
  95.     pm_error( "unusual bits per pixel (%d), can't read", h.bitpix );
  96.     }
  97.     if ( h.naxis != 2 && h.naxis != 3 )
  98.     pm_error( "FITS file has %d axes, can't read", h.naxis );
  99.     cols = h.naxis1;
  100.     rows = h.naxis2;
  101.     if ( h.naxis == 2 )
  102.     images = 1;
  103.     else
  104.     images = h.naxis3;
  105.     if ( imagenum > images )
  106.     pm_error( "only %d image%s in this file",
  107.           images, images > 1 ? "s" : "" );
  108.     maxval = min( fmaxval, PGM_MAXMAXVAL );
  109.     if ( h.datamin != 0.0 || h.datamax != 1.0 )
  110.     scale = maxval / ( h.datamax - h.datamin );
  111.     else
  112.     {
  113.     /* FITS images are not required to contain DATAMIN and DATAMAX cards
  114.     ** In this case, it would be necessary to read through the image
  115.     ** twice to properly scale.  Take this shortcut instead, remembering
  116.     ** that all FITS integers are signed values.  This scheme works
  117.     ** on most images because most astronomical data reduction packages
  118.     ** scale the images when writing so as to make maximal use of the
  119.     ** dynamic range of the output format. */
  120.     h.bscale = 1.;
  121.     h.bzer = ( fmaxval + 1.0 ) * 0.5;
  122.     scale = maxval / fmaxval;
  123.     }
  124.  
  125.     pgm_writepgminit( stdout, cols, rows, maxval, 0 );
  126.     grayrow = pgm_allocrow( cols );
  127.     for ( image = 1; image <= imagenum; ++image )
  128.     {
  129.     if ( image != imagenum )
  130.         pm_message( "skipping image %d of %d", image, images );
  131.     else if ( images > 1 )
  132.         pm_message( "reading image %d of %d", image, images );
  133.     for ( row = 0; row < rows; ++row)
  134.         {
  135.         for ( col = 0, gP = grayrow; col < cols; ++col, ++gP )
  136.         {
  137.         int ich;
  138.         double val;
  139.  
  140.         switch ( h.bitpix )
  141.             {
  142.             case 8:
  143.             ich = getc( ifp );
  144.             if ( ich == EOF )
  145.             pm_error( "EOF / read error" );
  146.             val = ich;
  147.             break;
  148.  
  149.             case 16:
  150.             ich = getc( ifp );
  151.             if ( ich == EOF )
  152.             pm_error( "EOF / read error" );
  153.             val = ich << 8;
  154.             ich = getc( ifp );
  155.             if ( ich == EOF )
  156.             pm_error( "EOF / read error" );
  157.             val += ich;
  158.             break;
  159.  
  160.             case 32:
  161.             ich = getc( ifp );
  162.             if ( ich == EOF )
  163.             pm_error( "EOF / read error" );
  164.             val = ich << 24;
  165.             ich = getc( ifp );
  166.             if ( ich == EOF )
  167.             pm_error( "EOF / read error" );
  168.             val += ich << 16;
  169.             ich = getc( ifp );
  170.             if ( ich == EOF )
  171.             pm_error( "EOF / read error" );
  172.             val += ich << 8;
  173.             ich = getc( ifp );
  174.             if ( ich == EOF )
  175.             pm_error( "EOF / read error" );
  176.             val += ich;
  177.             break;
  178.  
  179.             default:
  180.             pm_error( "can't happen" );
  181.             }
  182.         *gP = (gray) ( scale * ( val * h.bscale + h.bzer - h.datamin) );
  183.         }
  184.         if ( image == imagenum )
  185.         pgm_writepgmrow( stdout, grayrow, cols, maxval, 0 );
  186.         }
  187.     }
  188.     pm_close( ifp );
  189.     pm_close( stdout );
  190.  
  191.     exit( 0 );
  192.     }
  193.  
  194. static void
  195. read_fits_header( fp, hP )
  196. FILE* fp;
  197. struct FITS_Header* hP;
  198.     {
  199.     char buf[80];
  200.     int seen_end;
  201.     int i;
  202.     char c;
  203.  
  204.     seen_end = 0;
  205.     hP->simple = 0;
  206.     hP->bzer = 0.0;
  207.     hP->bscale = 1.0;
  208.     hP->datamin = 0.0;
  209.     hP->datamax = 1.0;
  210.  
  211.     while ( ! seen_end )
  212.     for ( i = 0; i < 36; ++i )
  213.         {
  214.         read_card( fp, buf );
  215.  
  216.         if ( sscanf( buf, "SIMPLE = %c", &c ) == 1 )
  217.         {
  218.         if ( c == 'T' || c == 't' )
  219.             hP->simple = 1;
  220.         }
  221.         else if ( sscanf( buf, "BITPIX = %d", &(hP->bitpix) ) == 1 );
  222.         else if ( sscanf( buf, "NAXIS = %d", &(hP->naxis) ) == 1 );
  223.         else if ( sscanf( buf, "NAXIS1 = %d", &(hP->naxis1) ) == 1 );
  224.         else if ( sscanf( buf, "NAXIS2 = %d", &(hP->naxis2) ) == 1 );
  225.         else if ( sscanf( buf, "NAXIS3 = %d", &(hP->naxis3) ) == 1 );
  226.         else if ( sscanf( buf, "DATAMIN = %lf", &(hP->datamin) ) == 1 );
  227.         else if ( sscanf( buf, "DATAMAX = %lf", &(hP->datamax) ) == 1 );
  228.         else if ( sscanf( buf, "BZERO = %lf", &(hP->bzer) ) == 1 );
  229.         else if ( sscanf( buf, "BSCALE = %lf", &(hP->bscale) ) == 1 );
  230.         else if ( strncmp( buf, "END ", 4 ) == 0 ) seen_end = 1;
  231.         }
  232.     }
  233.  
  234. static void
  235. read_card( fp, buf )
  236. FILE* fp;
  237. char* buf;
  238.     {
  239.     if ( fread( buf, 1, 80, fp ) == 0 )
  240.     pm_error( "error reading header" );
  241.     }
  242.